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Abstract. We consider a reduced order model of an air-breathing hypersonic engine with a time- 
dependent stochastic inflow that may cause the failure of the engine. The probability of failure is 
analyzed by the Freidlin-Wentzell theory, the large deviation principle for finite dimensional stochastic 
differential equations. We compute the asymptotic failure probability by numerically solving the 
constrained optimization related to the large deviation problem. A large-deviation-based importance 
sampling suggested by the most probable inflow perturbation is also implemented to compute the 
probability of failure of the engine. The numerical simulations show that the importance sampling 
method is much more efficient than the basic Monte Carlo method. 
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1. Introduction. Operation of a scramjet (supersonic combustion ramjet) is 
difficult to accurately model with state-of-the-art codes due to the complex physical 
and chemical systems governing the flow of air through and the combustion in the 
engine. When coupled with environmental uncertainties, a priori design of a safety, 
high-performance operation plan (fuel schedule) is a lofty ambition. This paper nu- 
merically analyzes the probability of failure for a scaled engine operating at about 
Mach 2. While there are many uncertainties outside of the engine that can effect 
the operability, we focus only on those within the core engine-system: the isolator, 
combustor and nozzle. The combustor is the location of the most complex (and least 
certain) chemistry: the amount of heat released here directly effects the performance 
of the system. 

Thermal choking can result from excessive fueling, which decreases immediate 
performance and produces a shock that may travel through the isolator. If the shock 
reaches the entrance of the isolator, the engine will stall; this is called unstart. While 
there are many other causes of unstart (for example, thermal deformation of the engine 
PQ), we focus on the unstart directly related to the inflow perturbations [TT] and the 
fueling of the engine. Iaccarino et al [5] studied the relationship between the amount 
of heat released and the operability of the engine over short time scales and with 
a low-fidelity model of the stochastic nature of the fueling. This paper studies the 
resulting uncertainty of the flow in the engine with a stochastic inflow Mach number. 

The Euler equations are frequently employed to model compressible flows in 
aerospace models, especially in steady-state as a numerically tractable model when 
designing airfoils (see [5]). The time-dependent, quasi-one-dimensional form of the 
equations can be used to capture the geometry of a compression-expansion engine 
(see e.g. [H]) and replicates many of the physical features of actual flows in scramjet 
engines. When used to model scramjets, a forcing term is added that models the heat 
release mechanism of fueling. It is well documented that these models capture the 
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physical shock that results from thermal choking. Due to the one-dimensional nature 
of these equations, they can be solved quickly and are ideal as a reduced order model. 

The large deviation principle is used to analyze events with exponentially small 
probability. The Freidlin-Wentzell theory, the large deviation principle for finite di- 
mensional stochastic differential equations is the mathematical tool to compute the 
probability of failure: when the random perturbation in the stochastic differential 
equation is small, the probability of failure decreases exponentially fast and the rate 
of decay of the probability is governed by the minimum of the rate function over the 
event of interest. We numerically solve this constrained optimization problem to ob- 
tain the asymptotic probability of failure and the most probable path causing unstart 
under several interesting cases. This so-called the minimum action method has been 
successfully applied to different model problems (see [4l 116]). 

Another obvious way to compute the probability of unstart is to use the Monte 
Carlo simulations. It has been extensively used in the engineering community to con- 
sider more elaborate scramjet models such as the two-dimensional model with the 
second order discretization, and it can be accelerated by using the adjoint-based sam- 
pling method [14] . When the targeted probability is small, however, because of the 
natural limitation of the basic Monte Carlo method, one needs a excessively large 
number of samples to accurately estimate the probability, and such large amount of 
computations leads to the inefficiency of the basic Monte Carlo method. We use 
the large-deviation-based importance sampling technique, the importance sampling 
Monte Carlo method whose change of measure based on the minimizer of the large 
deviation principle. Our numerical results show that the large-deviation-based im- 
portance sampling outperforms the basic Monte Carlo method. 

This paper is organized as follows: in Section [2] we discuss the equations that 
govern the flow, the geometry of the engine and the definition of the unstart; Section 
[3] briefly introduces the classical Freidlin-Wentzell theory, the large deviation principle 
for finite dimensional stochastic differential equations, and the large deviations for the 
related Euler schemes. In Section [2J we explain how to formulate the unstart of the 
scramjet as a large deviation problem. Section [5] shows the numerical results of the 
large deviation problems in Section [4] under different settings. In Section[6]we use the 
importance sampling Monte Carlo method based on the solution of the large deviation 
problems in Section [5] to directly estimate the probability of the unstart. Section [7] 
concludes this paper. The table of parameters and the numerical PDE method for 
the governing equation are in the appendices. 

2. Model Problem. 

2.1. Governing Equations and Engine Geometry. The quasi-lD compress- 
ible Euler equations serve as our reduced model of the engine-combustion system and 
capture the phenomena of unstart due to fueling. This model was developed by Iac- 
carino et al [5j and is similar to the model developed by Bussing and Murmam [3J. 
The quasi-lD compressible Euler equations are the following hyperbolic system: 

pu 

pu 2 + P I ='-^R\\P\ -\pu 2 + P\\ + \ I, (2.1) 



i 



K (E + P)uj 




where p is the density of the fluid, pu is the momentum and E is the total energy. The 
pressure P can be derived from an equation of state and we take P = (7— 1) (E— pu 2 /2) 
where 7 is the ratio of specific heats, taken to be 1.4. The function A(x) describes the 
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Fig. 2.1. Geometry of ID engine model. This figure is a schematic of both the engine 
geometry (shown in the lower half of the image) and of the fueling profile and a resulting shock. 
Lj, Lq and Lg are the lengths of the isolator, combustor and expansion region; 9j, 9q and (9g are 
their respective angles; the air flows in from the left and out at the right. In the top half, the x axis 
gives the location in the engine and the y axis is time. T is the full length of a fueling period and b 
is the burst length. 



cross-sectional area of the engine; we assume that the width is constant and thus the 
area varies as the height; see Figure [2TT] for a sample height profile and the following 
mathematical definition: 

!A — xsinOj, —Li < x < 0, 

A Q + xsm9 c , 0<x<L c , 
A + L c sin 6 C + (x - L c ) sin 6 E , L c < x < L c + L E - 

The term f(x,t) models the heat release due to fueling and takes the form: 



f(x t) - H^f^ ' ^ ' f stoch ' H P r °P ' A ° ' P ° ' U °/( L C ' A ( x )): x € [0,Lc], 



I 



otherwise. 



(2.2) 



where f(x) = x 1 / 3 following Iaccarino et al [5] and Riggins et al [IU] and O'Byrne 
et al jS] . <t> is the equivalence ratio and governs the amount of heat released into the 
system per unit time and f(t) is an indicator for when the engine is fueling. 



2.2. Unstart of the Engine. The Mach number, defined as M = u/^^/P/p, 
characterizes the behavior of the flow. The Mach number of the flow is plotted Figures 
2.2(a)|2~2(d)| for different values of 0, when the engine is fueled from t = 0.5ms to 



t = 1.5ms. In all cases a shock develops, from a supersonic Mach number (Mach 2, 
green) to a subsonic Mach number (0.4, blue). This shock extends into the isolator of 
the engine and persists after the fueling has stopped. The distance into the isolator 
that the shock travels before receding and the amount of time it takes for the engine 
to return to a "normal" idle state are functions of how much head is injected. In 



Figure 2.2(d) the shock reaches the left boundary of the inlet; this is called "unstart" 
and the engine has failed and ceases to produce thrust. To determine if the engine 
has unstarted, the shock location, defined as: 



Xshock = sup{a; G [-L c , 0] : M(x) > 1}, 
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Fig. 2.2. Effect of heat release parameters on deterministic solutions of Mach number, (a)- 
(c):Inflow Mach number = 2.0, b = lms. Progressing from (a) to (c), more head is added per unit 
time, driving the shock further into the isolator. In (d), where the shock reaches the left boundary, 
the engine stalls; this is called unstart. (d) Here, T/b = 3.33; the bursts are not spaced sufficiently 
far apart, and the location of the shock compounds. 



is tracked. When the engine is in a idle state x 



shock 



0; when the engine has failed 



shock 



2.3. Heat Release Model. With this simplified model, there are two regimes 
of operation depending on 4>: for sufficiently small (f> (Figure 2.2(a) I the shock does 
not leave the combustor, however insufficient thrust is produced; for greater than 
a threshold value <jj* w 0.25, a shock forms in the isolator and will propagate until 
the engine unstarts. This suggests that the simplest heat-release program that could 
result in sustained operation of the engine is to inject heat periodically. In this paper 
f(t) is defined as follows: 



/(*) 



f , t £ [nr, tit + b), 
0, te[riT + b, (n + l)r). 
4 



where r is the length of the fuel cycle and b is the length of the fuel burst (the 
amount of time fuel is being injected into the engine). With this heat release model, 
a potential cause of unstart is not spacing the fuel burst sufficiently far apart (t jb too 
small) which causes the shock to build upon itself and eventually leads to unstart, see 
Figure |2. 2(d)] 



The instantaneous thrust produced by the engine is given by 

thrust(i) = m e u e - m lUl + (P e - P t )A e = [Apu 2 ) e - (Apu 2 ), + (P e - Pi)A e , 

where the index e is for quantities at the exit of the engine and i if for quantities at 
the front of the inlet. The mass flow, m, is given by the state variable pu. When the 
engine is stalled no additional thrust is produced. The thrust produce is proportional 
to the amount of heat injected. 

Depending on the extensive numerical experiments in |15j , two fueling profiles are 
used in this paper: <fis = 0.78, t$ — 0.5ms, 65 = 0.1ms (the short fuel cycle) and 
4>l — 0.78, tl = 2ms, 6l = 0.4ms (the long fuel cycle). Note that because 4>s = 4>l 
and rs/bs — r^jb^ — 5, these two fueling profiles release the same amount of the 
heat so they generate roughly the same amount of thrust. 

2.4. Inflow Uncertainty and Its Effect on Unstart. There are many sources 
of uncertainty that contribute to the total uncertainty in the operability of the scram- 
jet engine. Those specific to the engine are the geometry (due to manufacturing errors 
or imperfections), the ratio of specific heats 7, the inflow conditions (the Mach num- 
ber, density and pressure of air entering the engine) and the combustion processes. 
In this paper we only address the uncertainty of the inflow Mach number Mi n (t) that 
is modeled as: 

M m (t) = M m (0) + ea M W t (2.3) 

with the initial condition Mi n (0) = 2, where e and a are positive constants and Wt is 
the standard Brownian motion. In addition, we assume that the inflow density pi n (t) 
and the inflow pressure Pi n (t) are constant in time. Then the perturbation of M in (t) 
completely comes from the perturbation of the inflow speed u in (t) 

The inflow uncertainty can greatly affect the stability of the screamjet even though 
it operas normally under steady inflows. As shown in Figure [2~3} if the aforementioned 
fueling profiles are used (the short and long fuel cycles) and M in (t) = M in (0) (ea = 0), 
the scramjet operates normally. However, if stochastic inflows are used (ea 7^ 0), then 
the inflow may affect the locations of the shocks and therefore the probability of the 
unstart is nonzero. 

In this paper, the primary goal is to consider the probability of the unstart due to 



the stochastic inflow modeled as (2.3). We are especially interested in the case that 



the parameter e in (2.3 1 is a small positive value, which lead to the large deviation 



analysis in the next section. 

3. Large Deviation Principle. In this section, we first briefly review the clas- 
sical Freidlin-Wentzell theory, the large deviation principle (LDP) for the stochastic 
differential equation. Then the analogous LDP for the Euler scheme of the discrete 
problem is introduced. The discrete problem provides a good understanding of the 
LDP for our unstart problem in the next section. 

3.1. The Freidlin-Wentzell Theory. We consider the following stochastic dif- 
ferential equation (SDE): 



dX e t = a{Xl)dt + eadWt, 
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(3.1) 
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Fig. 2.3. Effects of the inflow uncertainty on the scramjet stability. The left figures are the Mach 
numbers in the engine with constant inflow (Mach 2). The engine operates normally under both 
fueling profiles (the short and long cycles) because the subsonic flows do not reach the left boundary. 
After replacing the steady inflows by stochastic ones (shown in the right figures) , however, the shock 
locations are affected by the inflow and reach the left boundary so the unstart happens (shown in the 
middle figures). 



where Xq = xq is deterministic, Wt is the standard Brownian motion, and e and a 
are positive constants. The standard theorem says that there exists a unique strong 
solution on any finite time interval [0, T] if a : R — > K is uniformly Lipschitz continuous 
(see, for example, [5]). We assume that X t on [0,T] with Xq — xq is the unique 
solution of the following differential equation: 

jX t = a(X t ). 

It is well-known that as e — ¥ 0, X£ converges to X t in probability, that is, for any 
5 > 0, 



lim P max \Xf - X t \> 5) = 0. 

e^O \te[0,T]' ) 

The Freidlin-Wentzell theory says that as e — > 0, the asymptotic probability that 
X\ deviates from X t can be computed as follows: 

- inf l(x t ) < liminf \ logP(X f e 6 A) 

< lim sup \ logP(X t e EA)<- inf X(x t ), 

where the rare event A is any subset of C Xo ([0,T]), the space of continuous paths 
on [0, T] with the starting point xq endowed with the standard topology. A and A 
stand for the interior and closure of A, respectively. The rate function X measures 
how much X t deviates from its typical behavior in the L 2 sense: 

1{xt)= \£*fi[it x t-<x(xtj?dt, x t =x + J t y(s)ds, yeL 2 ([0,T}), 
I +oo, otherwise. 



Remark. For the proof of the Freidlin-Wentzell theory, readers may refer to [3| 
Section 5.6]. In fact, the complete version of the Freidlin-Wentzell theory works for 
any finite dimensional SDEs and a can be a function of X t . However, because we 
simply model the inflow Mach number as a Brownian motion, we use the current 
version to keep the expression of the rate function simple. 

An immediate observation is that if A is regular, that is, 

inf T{xt) = inf I[xt) = inf I(xt), 
x t eA x t^ A x t eA 

then we can represent the asymptotic probability in the following way: 

P(A7 e A) ~ exp f -\ inf X(xt) 

for small e. In addition, if X t £ A, then P(X t e £ A) ~ 1 as e — >• 0. Indeed, because 
XI — > X t as e — > 0, any regular event containing X t should be of probability one as 
e 0. 

3.2. Large Deviations for the Euler Scheme. For the computational pur- 
pose, one can use the Freidlin-Wentzell theory to obtain the analytical rate function, 
then discretizes the rate function and the functional space, and finally obtains the 
probability by solving the numerical optimization. However, it is more informative 
to consider the discrete problem by the Euler method and derive the exact LDP 
for the discrete problem. Then this LDP can be solved directly by the numerical 
optimization. 

We consider the following difference equation by the Euler scheme: 

K+i = K + + eaAW n+1 , (3.2) 

where Xq = xq is deterministic, {AW n +i}n=o are independent Gaussian random 
variables with mean zero and variance At, and e > 0. We assume that T — NAt. 
The joint density p of {X{, . . . , Xf^) can be derived by the Markov property: 

AT-l 

P( X N =X N ,...,Xl= Xi) = Y[ h i X n+l = x n+l\X e n = X n ) 

n=0 

N ~ 1 1 / 1 \ 

= S exp ~ Xn ~ a ^ At f) 

= (2^Atr^ exp (^±^ - a(x n ) 

\ n=0 ^ 

Let x = (xi, . . . , xpf) and 




Then given a set A c R N the probability that (X{, . . . , X e N ) e A is 

P(K,.,4)£A) = j (2Tre 2 a 2 At)- N / 2 exp(-^I(x)^dx. (3.3) 
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We use Laplace's method to compute the asymptotic probability as e — > 0. The 
basic idea is that as e — > 0, the mass of the integrand in (3.3) will concentrate at 
its maximizer, which is the minimizer of /. We can therefore use the minimum to 
compute the probability. 

Theorem 3.1. (Laplace's method) Assume that x* — arg min^g a I (%) with x* G 
A, and I(x) grows at least quadratically. Then 

lime 2 log P^,...,^) G A) =-I(x*). 

Proof. Because I(x) grows at least quadratically, for any 5 > 0, we can find a 
sufficiently large ball Bd = {x G : ||x|| < d} such that 



J {2ire 2 At)- N / 2 exp(-^I{x)^dx<5 



for all sufficiently small e. To compute the upper bound, it hence suffices to prove the 
case that A is bounded. We rewrite the integral as 

exp (-~I{x*)\ J (2^ 2 A<)- JV / 2 exp f-^[I{x) - I(x*)]\ dx. 

The upper bound is 

e 2 \ogP((Xt,...,X%) G A) 

= -I(x*) +e 2 log{2ne 2 At)- N / 2 + e 2 log j exp (-^{I(x) - I(x*)]j dx 

< -I(x*) + e 2 log(2 7 re 2 Af)- Ar / 2 + e 2 log |A|. 
Then we have the upper bound of the limit 

lim e 2 logP((X 1 £ , . . . , Xff) G A) < -I(x*). 

To show the lower bound, we use the assumption that x* G A. By the continuity 
of /, for any 5 > 0, there exists a neighborhood N$ c A of x* such that < 
I(x) — I(x*) < S for x G Ng. Then we have the lower bound: 

e 2 log¥((Xf, ...,4)eA)> e 2 log¥((Xt, .. .,X* N ) G N s ) 

= -I(x*)+e 2 \og(2TTe 2 At)- N ^ 2 + e 2 log J exp (~^[I(x) - I{x*)]j dx 

> -I(x*) + e 2 \og(2ne 2 At)- N / 2 + e 2 log(-\N 5 d). 
Therefore the lower bound is also obtained: 

lim e 2 logP((Xf, . . . , X%) G A) > -J(z*)- 

□ 

Readers can find that I{x) — > X{xt) as At — > 0. In fact, the following lemma says 
that X* is an exponentially good approximation of X\ . 

Lemma 3.2. (Exponential equivalence \S, Lemma 5.6.9]) If W t in (3.1) and 
AW n+ i in (Oj satisfy 



An+l)At 

AW n+1 = dW s = W (n+1)At - W nAt 

JnAt 



almost surely, then for any S > 0, 




lim limsuplogP max sup |X t e — > S 1 = — oo. 

Af^O £ ^ V n te[nAf,(n+l)At] " / 

Lemma [372] also tells us that we can effectively simulate a rare event by the Euler 
method. The standard theory shows that the Euler method has the order of accuracy 
y/Ai; however, by the large deviations, the probability of a rare event is of order 
exp(— 1/e 2 ) for small e. Therefore we might need to decrease At exponentially in e 
to obtain the correct simulation, and such exponential discretization will make the 
simulation computationally impossible. Thanks to Lemma [372} we only need to choose 
a uniform At and the Euler method can accurately simulate a rare event for all small 
e. 

4. Large Deviations for the Unstart. In this section, we derive the large 
deviation principle for the unstart of the scramjet. First we formulate the analytical 
problem. Readers can immediately identify that the analytical problem is an appli- 
cation of the Freidlin-Wentzell theory. However, as this problem can only be solved 
numerically, we also need the LDP for the discrctized problem. In other words, we 
derive the LDP for the numerical PDE of the flow equation. The readers can also 
find that the theory is almost ready because of the derivation in the last section. 

4.1. Analytical Model. Recall that the flow equation in the scramjet engine is 
p\ ( Pu 

pu\ + pu 2 + P | = ^hf ( ( P | ( pu 2 + P | | + ( | , (4.1) 
Ej t \(E + P)uj_ 

for x £ [— Lj, Lc + Le] and t £ [0, T]. To solve this PDE, one also needs to specify 
the initial condition and the inflow boundary condition. We assume that a suitable, 
deterministic initial condition is given. For simplicity, suppose that P(t, —L{] = Pq 
and pit, —Li) = po for all t € [0, T]; then the perturbation of the inflow Mach number 
M in (t) is entirely from that of the inflow speed u,* n (t) := u(t, —Lj) as M = u/ y/^jPjp. 
The inflow speed is modeled as a Brownian motion: 

Ui„(t) = u m (0) + ea u W t , 

with Uj„(0) = uo = 1300(m/s) that corresponds to Mach 2 in our setting. Then from 
Section [3~ij Ui n {t) satisfies the large deviation principle with the rate function 

T( A- i 2^ Jo (ft U *n(t)) 2 dt, Mln W=«m(0)+/„'^,l/eI 2 ([0,T]), 
-L(Uin) — \ 

l+oo, otherwise. 

by the Freidlin-Wentzell theory. The rare event A is the set of all possible Ui n (t) that 
causes the unstart in the time horizon [0, T], the event that the subsonic flow reaches 
the entrance of the isolator during [0, T]: 

A = {u in {t) : u in (0) = u , 3t € [0,T], x shock (t) = -£/}. 

Then the probability of the unstart is obtained by the large deviation principle: 

F(u in £ A) ~ exp ( -\ inf X(u in ) 



for small e. We note that in the optimization problem inf Uin ^X(ui n ), although the 
objective function T{u in ) is convex, it is very difficult to verify the convexity of the 
constraint set A, because it involves the analysis of the nonlinear hyperbolic system 



(4.1) 



4.2. Numerical Model. 

4.2.1. Numerical PDE. We uniformly discretize the space and time: —Lj = 
xq < ■ ■ ■ < xk — Lc + Le and = to < ■ ■ ■ < ijv = T. As the rate function needs 
to be defined a priori, we choose a uniform At satisfying the CFL condition so that 
the numerical PDE method is stable when we solve the large deviation problem. By 
convention, XJ! denotes the average of the quantity X over the cell (xk, Xk+i) at time 
t n . The local-Lax- Fricdrichs (LLF) scheme is used to solve numerically the governing 



equation (2.1 1. See also Appendix B for the details of the numerical PDE method. 



4.2.2. The Rate Function. In the discrete case, similarly, we model the inflow 
speed Uinin) as a Gaussian random walk: 

u m (n + l) = u in {n) + ea u AW n+ i, n = 0, . . . , N - 1 (4.2) 

with Ui n (0) = Uq = 1300(m/s). {AWn+i}^® are independent Gaussian random 



variables with mean zero and variance At. From Section 3.2 Uj„(n) satisfies the large 
deviation principle with the rate function: 



I{Ui n ) 



At ( u m (n + 1) - u in (n)' s ' 



2a?, ^ \ At 

u n=0 



We note that I(ui n ) is a convex function in it ira . 

4.2.3. The Rare Event. For the discrete problem, we define the unstart as the 
event that the subsonic flow is produced at (x\, X2), the location next to the entrance 
of the isolator (xq, x\). Hence the rare event A is the set of Uj„ = (uq, . . . , Uq ) causing 
the unstart on [0, T]: 

A = {u ln - (u° , ...,u Y ): u° = uq, min M™ < 1}, (4.3) 

l<n<N 

where M™ is the Mach number on (x\, x%) at time t n and is computed by the numerical 
PDE method. Similar to the analytical case, although the rate function is convex, it 
is difficult to verify the convexity of A. 

4.3. Numerical Optimization. By the large deviation principle, the probabil- 
ity of the unstart is therefore 

P(w m £ A)~ exp ( — \ inf I{u m ) 
\ e Mine a 

for small e. We compute this probability by solving the nonlinear constrained opti- 
mization problem inf„ inS A I(ui n ). We use the interior-point method (see Chapter 
19]) as the optimization algorithm. 

Here we address an important issue about this optimization. Any gradient-based 
optimization algorithm, including the interior-point method, only obtains a local min- 
imum, which might not be the global minimum. The global minimum is guaranteed 
only if the given problem is convex: both the objective function and the constraint 
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set are convex. In this problem, although the objective function I(ui n ) is convex, it 
is not clear if the constraint set A is. Thus in theory, we do not know if the answer 
we obtain is truly global. However, we solved this problem by multiple random initial 
guesses and obtained essentially the same result, therefore we have a high confidence 
that our answer is the global minimum. 

4.3.1. Dimension Reduction of the Optimization. Another computational 
challenge of the optimization problem is its extremely large scale. The typical setting 
of this problem is as follows. The space domain [— Lj, Lc + Le] is discretized into 100 
points. As the engine operates at the supersonic speed around Mach 2, the CFL con- 
dition requires that At s» 10~ 6 seconds. If we simulate the model up to 0.01 seconds, 
then we need 10 4 points to discretize the time domain. Every time we evaluate the 
constraint, a complete run of the numerical PDE is performed, and therefore every 
iteration of the optimization algorithm is very expensive. Consequently, to optimize 
the rate function with 10 4 variables is almost computationally impossible. 

To speed up the optimization, we can reduce the degree of freedom of the inflow 
speed Ui„. Indeed, although the CFL condition asks for the very fine grids in time, the 
resolution of the large deviation solution can be far lower than that. More precisely, 
we let N £ N such that mN — N for some m £ N, and Ui n (G), Ui n (m), . . . , Ui n {N) are 
the Gaussian random walks: 

Uin((n+l)m) = u in (nm) + ea u AW n+ i, n = Q,...,N-l, (4.4) 

where {AW n +i}n=o are independent Gaussian random variables with mean zero and 
variance mAt. The corresponding rate function is 

mAt ^-^ fu in ((n + l)m) - u in {nm)\ 2 
I(u in ) = -jjr L { J • ( 4 - 5 ) 

" n=0 V ' 

When we perform the numerical PDE, we linearly interpolate the other variables to 
obtain the inflow condition: 

/ k \ k 

u in {nm + k) = 1 I u in (nm) H u in ((n + l)m), 1 < k < m — 1. (4-6) 

\ m J m 

The optimization cost is generally at least proportional to the number of variables. 
In our case, we choose N = 20 and linearly interpolate the others, and therefore the 
algorithm only has to optimize over 20 variables. The reduced problem is at least 500 
time faster than the full problem. A potential issue of the reduced problem is that 
if the result of the reduced problem is accurate enough. Our numerical experiment 
shows that N = 20 is sufficient for this problem. In fact, in the later section, we find 
that even thought we double the resolution (N = 40), the relative improvement is less 
than 1%. 

Remark. By using the dimension reduction technique in this subsection, it is 
also possible to use the typical adaptive-time-discretization to solve the numerical 



PDE (4.1) by considering the following inflow condition: 



Ui n (t) = (l ] Ui n (nm) H Ui n ((n + l)m), t £ (nmAt, (n + l)mAt). 

\ m J m 

However, our numerical experiments show that the constraint set A in this way is less 
smooth that of the uniform discretization and the resulting numerical optimization is 
also less robust. For this reason, we still use the uniform discretization and choose an 
appropriate At. 
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4.3.2. Summary of the Numerical Optimization. We finish this section by 
summarizing the numerical optimization. 
Algorithm 4.1. 

1. The goal of this section is to compute numerically the optimization problem 
inf« ln eA^(um)- 

2. The full problem is computationally expensive so we instead solve the reduced 
problem 'm£u in eA I{u in ). 

3. The interior-point method is the optimization algorithm. The objective func- 



tion is given in (4-5) and the numerical PDE of (4-D is used to determine if 

^in £ A . 

4- Although our numerical experiments show that the optimization algorithm 
finds the same result even with different random initial guesses, a good initial 
guess can greatly speed up the optimization process. We let the initial guess 
{ui n (n)}^ =Q linear in n and choose Ui n (N) to be the largest value so that 
{uin(n)}n = o triggers the unstart. 

5. Numerical Results of the Large Deviations. In this section, we show 
the numerical results of the large deviation problem inf^.^gA I(uin)- The values of 
the parameters are listed in Appendix [A"| 

5.1. The Event of Subsonic Inflows. Before we go to the details of the large 
deviation result. We address an important subset of A that causes the unstart: the 
situation that the scramjet receives a subsonic inflow. We define 

B={u in (t); u in (0)=u o , min M m (t) = mm u in (t)/^/jP /po < !}■ (5-1) 

te[o,r] *e[o,T] 

We note that B C A and thus ml Uin ^X(ui n ) < inf Ui7iS B T(ui n ). Because itj n (i) is 
modeled as a Brownian motion, inf u . ne e^( u m) can be solved explicitly: 

u *n = arg Jnf^Cuin) = M - u + ^y/jPo/po = (l- ^ju + ^u . (5.2) 

In other words, the minimizer u* n is a straight line starting from uo and ending at 
0.5i*o, and equivalently, the minimizer M* n is a straight line starting from 2 and ending 
at 1. Further, we obtain an upper bound for inf Ui7ie ^I(ui n ): 

inf l{u in ) < inf l{u in ) = l(u* n ) = 7 — -u* n (t) dt = (5.3) 

Thus we can conclude that when inf^.^gA I{ui n ) is very close to this upper bound, 
the scramjet operates in a very stable region, because the probability of the unstart 
is as low as the theoretical limit. 

Remark. One may analogously want to define B in the discrete sense. However, 



in general B is not necessarily a subset of A because in the definition of A (see < 4.3 1), 
what we check is if the Mach number in the second cell MJ 1 is less than 1, but the 
inflow Mach number Mi n (n) ~ Mq < 1 does not imply M™ < 1. In fact one often 
needs Mi n (n) — Mq < 1 — rj with a small positive 77 to have M[ L < 1. Therefore in the 
later numerical results, readers will find that the computed infs in gA I{ui n ) in some 
cases will slightly higher than the continuum upper bound Uq/(8<t^T). Obviously this 
discrepancy will be reduced as we refine the discretization, and we still use this upper 



bound (5.3) as a reference value in our numerical results. 
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Fig. 5.1. The Mach numbers in the scramjet engine when a constant inflow (Mach 2) is used. 
The engine operates normally for both short and long fuel cycles. 





Short Fuel Cycle 


Long Fuel Cycle 






0.21504 


0.15603 


0.21125 



Table 5.1 

The optimal values of the rate function for the short and long fuel cycles. 



5.2. Impact of Fuel Cycles. The function of the fueling control is 

f(t,x)=c-<f>- /(*)• f(x), 

where c is a constant and f(x) is a fixed spatial function. <f> is the mixing ratio and 
f(t) is the control of the fuel cycle: 



/(*) 



1, te [ir, it + b], i € N, 
0, otherwise. 



In every fuel cycle of length r seconds, the fuel is injected in the first b seconds. 

Two representative fueling profiles are used: 4>s = 0.78, r§ = 0.5ms, b$ = 0.1ms 
and <f>L — 0.78, tj, = 2ms, b^ = 0.4ms. These two profiles are extensively studied 
in [15]. Because 4>s — 4>l and 6s/rs = 6l/tl, the two profiles generate roughly the 
same amount of thrust. However, the numerical experiments in |15j show that the 
first profile is more robust to prevent the engine from stalling. We observe the same 
qualitative behavior in the large deviation case. 



We can see in Figure 5.1 that with the constant inflow (Mach 2), the engine 
operates operates normally for both fuel cycles. No subsonic flows reach the entrance 
of the isolator. In this section, Figure [5TT] will serve as the reference case representing 
the normal situation. 

When the inflow perturbations are considered, we solve the large deviation prob- 



lem inf^. n6 A I(uin)- The solutions are plotted in Figure 5.2 (left). The optimal / 
of the short fuel cycle (0.21504) is about 1/3 higher than that of the long fuel cycle 
(0.15603). It means that when e is small, the probability of the unstart with the 
short fuel cycle case is lower than that with the long fueling case. This observation is 
qualitatively consistent with the study in |15) by the Monte Carlo simulations. 

Figure |5.2| also tells us an interesting information: with the short fuel cycle, the 
minimizer for inff liri£ A I{ui n ) is very close to a straight line starting from 2 and ending 
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Fig. 5.2. The large deviation solutions of two different fuel cycles. The optimal I of the short 
fuel cycle (4>s> T S> bs) is higher than that of the long fuel cycle (<j>L, bi,). It means the short 
fuel cycle is a more stable strategy. 



at 1, which is the minimizer u* n (t) for the large deviation problem inf Uin £sX(i 



discussed in Section 5.1 In addition, its optimal value of the rate function (0.21504) 
is also close to the continuum upper bound (0.21125), which says that with the short 
fuel cycle, the scramjet has very strong resistance to the unstart. 

5.3. Sensitivity to Constraints. We can slightly modify the constraint set A 
to see the sensitivity of the large deviation problem to the change of the constraint. 
More precisely, we define 



Arj.s 

Ai.2 



{Ui 



Uq 



uq, min M? < 0.8}, 

l<n<7V 

M , min M? < 1.2}. 

Kn<JV 



and let Ai.o := A in < 4.3 ) . It is reasonable to argue that Aq.s C Ai.q C A1.2 (by 



assuming that the Mach number is continuous in time.) Similarly, we can also define 
^0.8 j #1.0 and B1.2 according to (5.1), and derive the continuum upper bounds by the 



formulas similar to (5.2) and (5.3) (see also Table 5.2). 

From Figure |5.3| with the short fuel cycle, when we consider the constraint set 
A1.2, the result is not too surprising: the most probable inflow Mach number is close 
to a straight line starting from Mach 2 to Mach 1.2, and the shock generated by 
the engine has no contribution. Indeed, the set Ai .2 can be viewed as the case that 
the engine operates in a more normal condition, and in this case the effect of the 
engine is less than the effect in the case of A1.0 so the inflow Mach number is still 
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Fig. 5.3. The results of the large deviation problem by changing the constraint set Ai.n to 
l and Ai.2 with the short fuel cycle. 





Short Fuel Cycle 


Long Fuel Cycle 




inffi 4n eAo.8 I(v>in) 


0.26547 


0.18143 


0.3042 


inffit„eAi. 1{Uin) 


0.21504 


0.15603 


0.21125 


inffi in eAi. 2 1(u in ) 


0.13667 


0.13532 


0.1352 



Table 5.2 

The optimal rate functions over An. 8, Ai.o and A1.2. The values decrease with the looser 
constraints. 



the dominating effect. Nevertheless, when A0.8 is considered, the qualitative behavior 
changes. The most probable situation to have the Mach 0.8 subsonic flow is because 
the shock generated by the engine reaches the entrance, and the most probable inflow 
Mach number is not a straight line. In fact, the end point of the most probable 
inflow Mach number is higher than 0.8, which means that to in the sense of the large 
deviations, the scramjet does not need a inflow with Mach 0.8 to trigger the event of 
Ao.g. 

Then we consider the case of the long fuel cycle. From Figure |5~3} we find that for 
the case of A1.2, the most probable inflow Mach number (upper right) is close to that 
of the short fuel cycle case (upper right in Figure 5.3), and their optimal values of 
the rate function are also close to the continuum upper bound (see Table 5.2). That 
means the shock generated by the engine is not the dominating effect in this case. For 
the case of An,. 8, we can see that to have a subsonic flow of Mach 0.8 at the entrance 
one just needs a supersonic inflow. 
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Fig. 5.4. The results of the large deviation problem by changing the constraint set 
i and Ai.2 with the long fuel cycle. 



Ai.n to 



From Table |5.2[ we can find that the optimal values of the rate function with the 
short fuel cycle are consistently higher than those with the long fuel cycle. This tells 
us that the strategy of the short fuel cycle is a more robust profile for the safety of 
the operation of the scramjet. 

5.4. Impact of Engine Geometry. As Iaccarino et al. indicates in [5j, the 
engine geometry greatly influences the safe operating region against the unstart. In 
theory, for a larger angle of combustor 9c, the engine can accommodate more heat and 
therefore the scramjet has a larger safe operation region. However, [5] also mentions 
that a excessively large 9c will leads to the flow separation and cause the loss of 
thrust. The quasi-lD model cannot capture this phenomenon and therefore it can 
not be seen as well in our results. Readers may refer to [5] to see more details of the 
impact of 9 C . 

We also find the similar qualitative results in the large deviation sense. We 
change 0q from 7.5° (the typical setting in [S]) to two extreme values in [5], Be = 2.5° 
and 8c = 12°. In [5], the safe operation region of the scramjet is very small when 
9c = 2.5°, and 9c = 12° is considered as the largest angle of combustor without 
causing flow separation. 



From Figure 5.5 we see that when 9c = 2.5°, the scramjet is not able to contain 
as much heat as before. Then the shock is build up and reaches the entrance. The 
optimal values of the rate function are very low so the unstart can happen very easily 
when 9c — 2.5°. On the other hand from Figure [5T6| when 9c = 12°, the optimal 
values of the rate function of both fueling profiles are close to the continuum upper 
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Fig. 5.5. The large deviation solutions by changing 8^ from 7.5° to 2.5°. The optimal values 
of the rate function are significantly decreased. 



inf-tt iTlG Ai o I(^in) 


Short Fuel Cycle 


Long Fuel Cycle 




Q c = 2.5° 


0.088937 


0.046034 


0.21125 


e c = 7.5° 


0.21504 


0.15603 


0.21125 


C = 12° 


0.21505 


0.2147 


0.21125 



Table 5.3 

The optimal values of the rate function of different 8c ■ The higher 8q have the higher optimal 
values 



bound and the minimizer are similar to u* n . This is because for a larger 0c, the engine 
has more space for the generated heat and the shock is not easy to move upstream. 

5.5. Resolutions of Large Deviation Solutions. From the previous simula- 
tions, we are actually convinced that N = 20 is the sufficient resolution of the large 
deviation solution due to the smoothness and the strong linearity of the solutions. 
Here we support our claim that N = 20 is enough by doubling N to 40. If the results 
are very close, we can strongly believe our a rgument. 



We use the same settings in Section 5.2 but let N — 40 instead. From Figure 



5.7 



the optimal solutions are essentially the same, and the relative differences between 
the optimal values of the rate function are less than 1%. 

6. Monte Carlo Simulation with Importance Sampling. The large devia- 
tion results obtained in Section [5] is the exponential rate of decay of the probability, 
but are not the actual one. In this section we compute the probability of the unstart 
by using the Monte Carlo method. 
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Fig. 5.6. The large deviation solutions by changing 9^ from 7.5° to 12°. The optimal value 
of the high fuel cycle case is significantly increased. The optimal value of the short fuel cycle case 
is almost the same because the value already reaches the upper limit. 



i n f«,„eAi 


H^n) 


Short Fuel Cycle 


Long Fuel Cycle 




N = 


20 


0.21504 


0.15603 


0.21125 


N = 


40 


0.21503 


0.15587 


0.21125 



Table 5.4 

The optimal values of the rate function for different resolutions. The relative differences are 
less than 1%. 



6.1. Introduction to Importance Sampling. The most basic way to com- 
pute the probability of the unstart P(A) is to genera te m any i ndep endent sample 

(n)}n=o satisfy (jZil and Then use the 

to determine if vp. 



paths {u J in }j =1 , where for each j, { 
numerical PDE in Section 
estimator is 



4.2.1 



S A. The basic Monte Carlo 



pMC _ 



1 - 

3 = 1 



(6.1) 



Because E[P MC ] = P(A), P MC is an unbiased estimator, which means that P MC -> 
P(A) almost surely as J — > oo by the law of large numbers. In addition, by the central 
limit theorem, the error bar of P MC is proportional to its standard deviation 



Std(P A/c ) = ^j[ p ( A ) 
18 



2 (A)f/ 2 . 



Mach Number, Short Fuel Cycle 



Inflow Mach Number, Short Fuel Cycle 




-0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 

Space(m) 

Mach Number, Long Fuel Cycle 
High Resolution, min I = 0.15587 

<10 Alu 





0.004 0.006 

Time(s) 

Inflow Mach Number, Long Fuel Cycle 



-0.5 -0.4 -0.3 -0.2 -0.1 




Space(m) 



Time(s) 



Fig. 5.7. The large deviation solutions with the high resolution. The settings are the same as 



those in Figure 5.2 but the solution resolution N is doubled to 40. The results are nearly the same 



and the relative differences between the optimal values of the rate function are less than 1%. 



In order to have a meaningful estimate, the order of the error bar should not exceed 
the order of the estimated probability. Namely, 

Std(P MC ) J_ / _J_ _ \ 1/2 

P(A) VP(A) J [ ' 

The large deviation analysis tells us that as e < 1, P(A) decreases exponentially in 
e so at the same time J has to increase exponentially. The exponential growth of J 
will eventually cause the basic Monte Carlo method computationally impossible. 

To see what causes this computational difficulty, we note that the basic Monte 
Carlo estimator is simply the empirical frequency of Uj„ € A. When e is small, only 
a very small fraction of the samples is meaningful (in A), and most of the samples 
has no contribution to the estimation so resulting the inaccuracy of the estimator. 

The well-established method to solve this issue is to use the importance sampling 
technique. The idea is that since most of the samples under the original measure P 
have no contribution, we use a different measure Q to sample uj n so that there is 
a significant fraction of the samples contributing to the estimate. Since we bias the 
measure P, a correction is needed to obtain an unbiased estimate. More precisely, we 
have 

dP 

P(u m € A) = Ep[l A (Wt7i)] = E Q[WWm)^("m)], 
19 



where dP/dQ is the change of the measure. The importance sampling estimator is 

J 

■ I 



0=1 

where u\ n is sampled under Q. Note that P is unbiased as E[P IS ] = P(ui n G A) 
and its standard deviation is 

Std(P MC ) = 4=Std(l A (^„ Au m )). 

v J "<U! 

6.2. Large-Deviation-Based Importance Sampling. The main challenge of 
the importance sampling is how to choose Q (and therefore dP/dQ) to lower the 
standard deviation. A good choice may come from the solution of the large deviation 
problem inf^.^gA I(ui„). Assuming that the minimizer 

u*n = ar S inf A I(u m ), 

Uj„6A 

is unique, then for any open neighborhood TSS(u**) of u**, we have the following 
asymptotic conditional probability: 

V(u m G N(u**)\u in G A) = 1 - V(u ln G N c (uZ)\u m G A) 

= i - p(n c (^) n a)/p( A ) ■s 1 i - exp( "^; nf ;- gNC( ^ )f ; A/ f" )) n° i. 

exp(- ? mf fiineA -/ (Sin) ) 

In other words, the mass of the conditional probability is concentrated exponentially 
fast around the most probable path u**. This observation motivates us to choose the 
measure Q so that the sampled u\ 's are centered around u**. 



Now we construct Q and dV/dQ based on u**. Recall that under P, from (4.4 1 
we have 

n 

Uin((n + l)m) = u in (nm) + ea u AW n +i = u Q + ecr u 2J AWi+i, n = 0, . . . , N — 1, 
We let Q such that under Q, u in is a Gaussian random walk centered at u**: 

n 

u in ((n+l)m)=ut*((n+l)m) + ea u J2 A Wi+i, n = Q,...,N-l, (6.2) 

1=0 

where {AWn+ij^Tg 1 are independent Gaussian random variables with mean zero and 
variance mAt und er Q. The intermediate variables are still determined by the linear 
interpolation (4.6). Note that {ui n (nm)}^ =1 are jointly Gaussian under both P and 
Q so the change of measure can be obtained explicitly: 

dp exp (-j2^2-(u in - u ) T S" 1 (u i „ - u )J 

— (u in ;uZ) = t r-, (6.3) 

dQ exp (-^(u^ - u»)TS-i(u,„ - u-)) 

where Uo, u*^ and iii n are JV-dimcnsional column vectors: 

u = (u , . . .,u Q ), u** = (u**(m), . . .,u**(Nm)), u m = (u in (m), . . .,u in (Nm)), 
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and S is the covariance matrix of (AWi, . . . , AW^) under Q. 

In summary, the large-deviation-based importance sampling is implemented as 
follows: 

Algorithm 6.1. 



1. Compute it** = arg inf^eA I(v,i n ) nu meric ally by Algorithm 4-1 

2. Sample J independent u\ n under Q by (6.2). 



3. For each j, compute the change of measure ^(u J in ]u**) in (6.3) 
4- The large-deviation-based importance sampling estimator is 



6.3. Simulation Results. Here we test the two estimators of the probability 
of the unstart: the basic Monte Carlo estimator P MC and the large-deviation-based 
importance sampling estimator P IS . The event we test is P(tti„ G A) = P(ui n e Ai.o) 



(see (4.3) for the definition of A) with the short and long fuel cycles (see Section 5.2). 
As the inflow condition is random, instead of the uniform time increment At, we 
use the adaptive time increment to ensure that the CFL condition is satisfied (see 
Appendix [B] for more details.) 

We let J = 10 4 and test for e = 0.2,0.22,0.24, . . . ,0.4, where P(A) = 0(10^) 
for e = 0.4 and P(A) = O(10" 3 ) for e = 0.2. The (numerical) 99% confidence interval 
and the (numerical) relative error are two quantities measuring the performance of 
the estimators. We define the numerical standard deviations as 



>AfC\2 



J'=l 

(Stdf f = -l_^ ( l A( ^)^(^ ; ^)-P^) 



Then the numerical 99% confidence intervals are 

{ pMC _ ^g td MC pMC + ^ Std M C]) [P IS _ ^Stdf, P IS + ^Stdf], 

v J v J v J V J 

and the numerical relative errors are Stdf c /P MC and Stdj S /P /s . 



From Figure 6.1 and 6.2 we see that P IS has the more accurate estimates; the 
improvement of P MC increases when the estimated probability decreases. We also 
note that although the importance sampling is intrinsically designed for the small e 
situations, our simulations show that it also improves the non-small e cases. Because 
the rare of convergence of P MC and P IS are 1 / yj, the improvement of the speed is 
the square of the ratio of the standard deviations Stdj /C /Stdj' s . Then we see that 
the improvement of the speed of P 1 ranges from a factor of 4.2 (e = 0.4) to a factor 
of 381.5 (e = 0.2) when the short is used, and from a factor of 5.4 (e = 0.4) to a factor 
of 170.3 (e = 0.2) if we use the long fuel cycle. Roughly speaking, the ratio of the 
improvement of is proportional to 1/P(A). 

7. Conclusion. In this paper, we use the large deviation principle to analyze the 
probability of unstart of a scramjet due to the random perturbation of the inflow. The 
numerical analysis is performed under various comparisons: the impact of the fueling 
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Fig. 6.1. The basic Monte Carlo method and the importance sampling estimator for the 
probability of the unstart when the short fuel cycle is used (the top figures). The error bar are the 
99% confidence intervals. The bottom left figure is the plot of the estimated probabilities in the log 
scale. The bottom right figure is the plot of the relative errors Std* fC JpMC an( ^ stdj S /P IS , and 
the ratio o/ Stdj^ to Stdj 3 , which is the factor of the improvement of Stdj 5 . 

schedules, the sensitivity to the constraint sets and the effect of the engine geometry, 
some of which are also confirmed and are consistent in the previous literature using the 
Monte Carlo method (0 [T5]); namely, the central analysis and the large deviation 
analysis have the high consistency on the region of operation. Further, the large 
deviation analysis gives a sharper information by providing the most probable inflow 
perturbation that causes unstart. 

We also implement the large-deviation-based importance sampling to overcome 
the limitation of the basic Monte Carlo method. Our numerical results show that 
when the probability of unstart is small but still not negligible (for example, the 
order of 10 -3 or even 10~ 4 ), the importance sampling is significantly better than 
the basic Monte Carlo. The theoretical ratio of the improvement can be up to the 
reciprocal of the estimated probability. 
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Appendix A. List of Parameters. The listed parameters are the default 
values. We will mention the changes in the main text if different values are used. 



Variable 


Name 


Value 


Units 
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State Variable 








P 

u 

Uin 

pu 

E 
P 

M 

7 


f ^ ■ 

Density 

Flow Speed 

Inflow Speed 

Mass Flow 

Energy Density 

Pressure 

Mach Number 

Inflow Alach Number 

xxatio oi opeciuc xieato 


- 

1 A 

L .4: 


Kg j m 

m/s 

m/s 

kg/m 2 s 

J/m 3 

Pascals 


Geometry 








A 

Lj 

Lc 

Le 

Oi 

6c 

9 E 


Minimum Cross Sectional Area of Engine 

Inlet Length 

Combustor Length 

Expansion Region Length 

Angle of Inlet 

Angle of Combustor 

Angle of Expansion Region 


0.008 
0.5 
0.1 
0.1 
0.0 
7.5 
15.0 


m 2 
m 
m 
m 

Degrees 
Degrees 
Degrees 


Fueling 








9 
f 

J stoch 
-H-prop 

Po 

Uo 

Po 
rs 
bs 

Tl 

b L 


Mixing Ratio 

Stochiometirc Fuel /Air Ratio 
Fuel Heating Value 

Free Stream Density - Taken as Inflow Density 
Free Stream Velocity - Take as Inflow Velocity 
Free Stream Pressure - Take as Inflow Pressure 
Short Fuel Cycle Length 
Short Fuel Burst Length 
Long Fuel Cycle Length 
Long Fuel Burst Length 


C\ TO 

u.uzy 
1.2 x 10 8 
0.159 
1300.0 
47842.0 
0.5 
0.1 
2 
0.4 


J/Kg 

kg/m 3 

m/s 

Pascals 

ms 

ms 

ms 

ms 


Numerics 








T 

K 

N 

At 

N 


Terminal Time 
Number of Cells 

Number of Time Grids for Large Deviations 
Time Increment of the Euler schemes (3.2| ( 4.2 1 
Resolutions of Wi„ in (4.4| 
Volatility of Um 
Volatility of M in 


0.01 
100 
10 4 
lO" 6 
20 
10 4 
96.9020 


s 
s 

m/s 3 / 2 



Table A.l: Table of Numerical Values 



Appendix B. Numerical PDE Methods of the Governing Equation. In 



this section, we describe the numerical method to solve the governing equation (2.1 ) 




Given a spatial discretization: —h% = xq < ■ ■ ■ < xk — Lc + Le, X£ denotes the 
average of the quantity X over the cell (xk,Xk+i) at time t n . We use the component- 
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Fig. 6.2. The basic Monte Carlo method and the importance sampling estimator for the 
probability of the unstart when the long fuel cycle is used (the top figures). The error bar are the 
99% confidence intervals. The bottom left figure is the plot of the estimated probabilities in the log 
scale. The bottom right figure is the plot of the relative errors Std* fC JpMC an( ^ stdj S /P IS , and 
the ratio o/ Stdj^ to Stdj 3 , which is the factor of the improvement of Stdj 5 . 




wise, first order local Lax-Friedrichs (LLF) scheme [131 Algoritm 4.4] to solve (2.1): 
/ Pi \ h 

— Pk u 'k - T~ ("^+1/2 ~ Fk-1/2) 

V El j 

+ h A p±i4 ( ( p«) - ( P i{uif\ pz))+h( V (b.i) 

A{x k+1/2 ) yy Q j y iE n + P n H JJ \f( Xk+1/2l t n )J 

where x k+1/2 = 0.5(x fc + x k+1 ), I™ = (7 - - p£«) 2 /2). JT« ±1/2 are the 

numerical fluxes generated by the component-wise LLF method: 




- ^max{|c"+<|,|4 l + i 

where c k = U "fP k / p k is the speed of sound. Because the component-wise scheme 
works well for the low-order methods (see |13|). we use it to reduce the computational 
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cost. 

For the inflow conditions, we let p^ — po and Pq = Pq for all n, and Uq is 
governed by the stochastic inflow speed Uj„(n). For the outflow condition, the simple 
extrapolation is used: E% +1 ) = (pK-v u K-v E K-i)- 

The following strategy is used to find the suitable initial condition: we simu- 
late the numerical PDE with p(0,x) = p(t,—Lz) = po, it(0, x) = u(t,—Li) = uq, 
E(0, x) = E(t, —Li) = Eq and f(t, x) = 0. Then with the aforementioned outflow ex- 
trapolation, the numerical solution obtained by this setting has the equilibrium state 
(p e (x), u e (x), E e (x)) and use this state as the initial condition for the LDP. 

In Section [4] and [5] the time increment h is taken as the uniform constant At in 



(4.2). This is because the uniform time grid results in a smoother constraint set A 



in (4.3), which increase the robustness of the numerical optimization inf^ in GA I{ui n ). 
Because the inflow condition Ui n is a controlled variable, we can choose a sufficiently 
small At so that the numerical scheme is stable when we solve inf^ iiig A I{ui n )- On 
the other hand, in Section [6j as the inflow condition is random and not controlled, we 
use the adaptive time increment: 

Ji» = o.r Ax 



max fc |c£ + u£| 

to satisfy the CFL condition and avoid the instabilities due to extraordinary inflow 
conditions. 
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